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A  COMPARISON  OF  RADIATION  TRANSPORT  METHODS 
IN  AXISYMMETRIC  GEOMETRIES 


I.  INTRODUCTION 

Radiation  from  a  plasma  and  its  subsequent  transport  through  a 
surrounding  medium  is  a  fundamental  physical  process  in  several  experiments 
presently  ongoing  at  the  Naval  Reas'  ireh  Laboratory.  In  the  laser-target 
interaction  of  the  PHAROS  III  experiment,  radiation  plays  an  important  role 


in  heating 

the  backside  of  the 

target 

leading 

to 

the 

rearward  plasma 

blowoff . ^ 

Radiation  from  atomic 

lines 

is 

one 

of 

the 

few 

diagnostics 

available 

to  study  the  dynamics  of 

the 

hot , 

shocked 

cavity 

of 

ambient  gas 
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formed  by  the  frontside  blovoff.  One  aim  of  the  Z-pinch  of  the  GAMBLE  II 

facility  is  to  create  a  highly  energetic  burst  of  x-rays  at  maximum 

3 

compression  on  axis,  and  another  one  is  possibly  to  develop  an  x-ray  laser. 

Although  both  of  these  experiments  have  been  modeled  with  one 
dimensional  radiation-hydrodynamic  numerical  codes  (Refs,  above),  there  are 
certain  aspects  which  require  at  least  a  two  dimensional  treatment.  For 
instance,  the  plasma  blowoff  from  a  laser  target  is  asymmetric  about  the 
plane  of  the  target,  though  nearly  symmetric  about  the  laser  axis.  And  a 
sausage  or  Rayleigh-Taylor  instability  about  the  compression  axis  can  disrupt 
a  uniform  Z-pinch  which  is  required  for  maximum  efficiency. 

Furthermore,  both  experiments  involve  large  changes  in  length  scales 
from  the  initial  configuration  to  the  final  one.  In  the  laser-target 
interaction,  the  disturbed  plasma  expands  from  the  size  of  the  target,  a  few 
microns,  to  nearly  a  centimeter  at  blast  wave  formation.  In  the  Z-pinch,  the 
annular  gas  puff  starts  at  about  a  centimeter  in  radius  and  compresses  to 
less  than  a  tenth  of  this  radius.  Hence  a  Lagrangian  or  an  adaptive  spatial 
mesh  is  needed  to  follow  the  plasma  with  proper  resolution. 

This  report  addresses  the  problem  of  how  to  solve  the  radiation  transfer 
equation  in  an  axisymmetric  geometry  over  an  arbitrary  quadrilateral  mesh  in 

V1.inuM.ripi  .ipprovcil  Juno  ID.  1^X6 


the  R-Z  computational  plane.  Two  numerical  schemes  are  considered.  The 
first  scheme  is  a  generalization  of  the  Eddington  diffusion  approximation  for 
radiation  transfer  wherein  the  radiation  flux  is  proportional  to  the  gradient 
of  the  radiation  energy  density.  It  uses  a  simple  flux-limiter  to  constrain 
the  radiation  flux  to  physical  values  in  regions  of  steep  gradients.  The 
resulting  equation  for  the  radiation  energy  density  is  a  diffusive  one  and  is 
solved  by  a  matrix  inversion  technique  over  the  spatial  mesh.  The  second 
scheme  is  more  direct,  xt  solves  the  radiative  transfer  equation  along  a  set 
of  rays  in  three-dimensional  space  extending  from  each  cell  center  in  the 
distorted  two-dimensional  mesh.  The  resulting  intensities  are  then  averaged 
to  get  the  local  radiation  energy  density. 

The  present  ray  trace  scheme  is  a  generalization  of  the  one  proposed  by 

4  5  6 

Colombant  and  Winsor  ’  and  used  in  a  laser  target  model  .  Their  earlier 

approach  performed  a  crude  average  of  intensities  for  rays  outside  of  the 

computational  plane  and  was  restricted  to  square  meshes  in  that  plane. 

The  solutions  from  the  diffusion  and  ray  trace  schemes  are  compared 

against  a  closed  form  solution  of  the  radiative  transfer  equation  for  an 

emitting  cylindrical  core  in  Section  II.  The  schemes  are  also  compared 

against  each  other  for  problems  of  more  complex  configurations.  Each  scheme 

has  its  advantages  and  disadvantages  as  will  be  pointed  out  below  in  the 

Conclusions,  Section  III. 

Let  us  begin  by  discussing  two  different  formal  solutions  of  the 

radiative  transfer  equation.  Let  I(r,B, v,t)  be  the  specific  intensity  at 

position  r,  in  the  unit  vector  direction  Q,  at  frequency  v  and  time  t,  with 
2 

units  of  ergs/cm  • sec • Hz • s teradian .  Let  a  (r,v,t)  be  the  pure  absorption 

cl 

coefficient  with  units  of  cm  and  likewise  as(r,v,t)  the  pure  scattering 

coefficient.  Finally  let  e(r,\>,  t)  be  the  isotropic  emission  coefficient  with 

3 

units  of  ergs/cm  -sec-Hz.  Then  for  non-relativistic  flows,  the  transfer 
equation  is^ 


2 


g  ai<r^|','U  «  fl-7  I(t,S,v,t)  *  [oa(r,v,t)»5s(r,M,t)]I(r,fi,v,t) 


where  S*V  is  the  directional  derivative  along  fi,  and 


E(r,  v,  t)  =  i  J  dfi  I(r,  fi,  \>,  t)  (2) 

is  the  radiation  energy  density  per  unit  frequency. 

To  obtain  an  integral  formulatijn  of  the  transfer  equation  applicable  to 
any  coordinate  system,  consider  a  steady  state  situation  without  scattering. 
Then  rewrite  eqn.  (1)  through  a  transformation  of  variables  as 


-  ^  I(r-sQ,  Q,  m)  +  <ja(r,v)I(r-sO,Q,v)  .  , 


where  s  (  >  0  )  is  a  distance  backward  along  the  direction  fi.  Define  the 
optical  depth  as 


t(r, r-sfi, \>)  =  Jq  ds'  aa(r-s'Q,v)  ,  (3) 


and  use  the  integrating  factor  e  to  obtain  the  formal  solution 


I(r,C,  v)  -  I(r-sbfi,  v)  e 


-  -t(r,r-sb2, v) 


*  ,* b  T(r’r-S'B’V) 

J0  as  4rt  e 


,  (4) 


WE 


»■» 


where  s^  is  the  distance  from  the  point  of  observation  r  to  the  outer 
boundary  surface  along  the  direction  -2.  The  radiation  energy  density  at  r 
is  given  by  eqn.(2).  The  essence  of  the  ray  tracing  scheme  is  to  evaluate 
eqn.(4)  over  a  set  of  rays  {2}  for  each  observation  point  r. 

An  alternative  approach  which  includes  the  time  dependence  and 
scattering  is  to  take  angular  moments  of  the  transfer  equation  (1).  In 
cylindrical  coordinates,  dQ  =  sin9-d9-d4>, 


2  =  cos9  e  +  sin9,cos$  e  +  sin9‘sin$  e.  , 

z  r  ♦ 


and  the  transfer  equation  becomes' 


1  31  31  .  a  .31  ,, 

-  xr  +  sm9*cos$  +  cosQ 
c  9t  3r  T  3z 


r  I?  sin0‘sin*  *  < V<Ts)I 


4  ii  +  cas  4n 


In  Cartesian  geometries  Feautrier  variables  are  often  employed  to  transform 
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j mr'1  --  ’  But  in  an 


the  transfer  equation  into  a  simpler  diffusion  equation. 


axisymmetric  geometry  the  presence  of  the  angular  derivative  in  eqn.(6) 
negates  this  transformation.  Taking  the  zeroth  angular  moment  of  eqn.(6) 
gives 


If  +  *-F  +  caaE  =  c  ’ 


where  the  radiative  flux  is  defined  by 


F( r ,  v,  t )  =  J  d2  I ( r ,  B,  v,  t ) 2 
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It  is  well  known  that  the  moment  method  suffers  from  the  standard 


problem  of  non-closure,  for  note  that  the  next  angular  moment, 

Z  If  +  cy-P  +  <Vas)P  =  0  ’  (9) 

introduces  the  radiation  pressure  tensor 

P(r,v,t)  =  ^  /  dS  I(r, 8, v, t )  8  8  .  (10) 

The  essence  of  the  diffusion  scheme  is  to  cut  off  the  moment  development  by 
relating  P  -  fE,  v*>ere  f  is  the  so  called  tensor  Eddington  factor. 
Neglecting  the  time  derivative  in  eqn.(9)  shows  that  this  relation  leads  to 


F  = 


7*  f  E  . 


(11) 


_s 


i 


i 
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II.  COMPARISON  OF  NUMMERICAL  SOLUTIONS 


A.  Simple  Model  Test  Problem  and  Closed  Form  Solution. 

In  order  to  compare  the  diffusion  scheme  for  solving  the  radiative 

transfer  equation  with  the  ray  tracing  scheme  ve  consider  an  axisymmetric 

test  problem.  This  first  problem  should  be  simple  enough  that  an  accurate 

closed  form  solution  can  be  obtained,  yet  it  should  have  the  character  of  a 

two  dimensional  structure.  Let  the  computational  plane  within  which  the 

distorted  mesh  lies  extend  upward  1  cm  in  the  z-direction,  i.e.,  the  symmetry 

axis,  and  outward  1  cm  in  the  radial  direction.  Rotate  this  square  around  the 

symmetry  axis  to  obtain  the  computational  volume  depicted  in  Fig.l.  A 

3 

spatially  and  temporally  constant  emissivity  of  e  =  100  ergs/cm  ‘sec-Hz  is 
assumed  within  a  cylindrical  volume  of  height  1  cm  and  radial  extent  0.2  cm. 
The  absorption  coefficient  within  this  region  is  a  =  10.0  cm-*.  We  neglect 
scattering  processes.  Outside  of  this  region  there  is  no  source  (e  =  0)  and 
the  absorption  coefficient  is  0.1  cm-*.  Thus  the  center  of  the  emitting 

region  is  optically  thick,  while  the  source  free  region  is  purely  absorbing 
and  optically  thin.  Near  the  source  region  the  problem  is  nearly  planar,  but 
farther  out  in  radius,  and  near  the  edges  of  the  computational  volume,  the 
finite  extent  of  the  source  in  the  z-direction  has  its  effects.  We  seek  the 
steady  state  solution  to  this  problem  over  the  distorted  mesh  shown  in  Fig. 2. 

Due  to  the  simple  geometry  and  spatial  dependence  of  the  emissivity  and 
absoption  coefficients,  the  specific  intensity  at  an  arbitrary  point  r  in 
the  computational  plane  can  be  computed  exactly  from  eqn.(4)  along  a 
direction  Q.  The  vector  Q  is  defined  by  the  angles  9  and  $  of  Fig.l,  and 
from  trigonometric  relations  the  optical  depth  and  line  integral  of  eqn.(4) 


are  evaluated.  Then  to  obtain  the  radiation  energy  density  at  r,  the  solid 
angle  subtended  by  the  source  region  at  position  r  is  discretized  into  a 
large  number  of  9  and  $  angles.  The  same  procedure  is  done  for  the  radiative 
flux  of  eqn.(8).  Finally,  we  solve  for  the  angular  integrated  radiation 
quantities  at  the  center  of  each  cell  in  Fig. 2.  The  results  for  the 
radiation  energy  density  (ERAD)  and  the  normalized  z-component  of  the 
radiation  flux  (FZ/c-ERAD)  are  displayed  in  the  contour  plots  of  Figs. 3a  and 
3b,  respectively.  The  quantity  in  the  parentheses  at  the  top  of  Fig. 3a  is 
the  total  radiation  energy  in  the  whole  computational  volume. 

B.  Diffusive  Solution. 

The  essence  of  any  radiation  diffusion  scheme  is  a  good  estimate  of  the 
Eddington  factor  f  of  eqn.(ll).  One  way  to  estimate  f  is  to  use  the  results 
from  a  simple  ray  tracing  solution,  but  this  defeats  the  basic  purpose  of  a 
diffusion  approach  in  an  axisymmetric  geometry.  Several  derivations  have 
been  proposed  for  relating  F  directly  to  VE  with  a  coefficient  R  that  limits 
1 F | / c * E  below  the  streaming  limit  of  uni ty . ^ ^ ^ ^ ^  All  of  these 
derivations  take  the  specific  intensity  to  be  isotropic  about  the  direction 
of  the  radiative  flux,  and  neglect  the  spatial  and  temporal  dependence  of  the 
Eddington  factor  f.  Vith  these  assumptions  an  expression  for  R  can  be 
derived.  It  is  easy  to  see,  however,  that  in  axisymmetric  problems  the 
radiation  need  not  be  isotropic  about  the  flux  direction.  Furthermore,  as 
noted  by  the  above  cited  authors,  in  source  free  or  scattering  free  regions 
these  theories  predict  that  the  radiative  flux  immediately  attains  the 
streaming  limit,  i.e.,  they  break  down  in  a  purely  absorbing  medium.  The 
experiments  mentioned  in  the  Introduction  are  characterized  by  a  small 


scattering  coefficient  and  a  significant  region  of  minimal  emission  but 
strong  absorption.  Any  attempt  to  remedy  these  flux-limiting  theories  would 
be  ad  hoc,  so  we  use  the  simplest  flux  limiting  scheme,  which  itself  is  ad 
hoc,  namely, 


P  = 


■u  \  <>  I  TO  I 

3<Vas)+2  E 


VE 


(12) 


This  formula  was  used  in  early  versions  of  LASNEX.3^  In  the  limit  of  small 
gradients,  eqn.(12)  reduces  to  the  classical  Eddington  approximation.  In 
regions  of  steep  gradients  one  finds  |F|  =  cE/2,  the  proper  result  for  the 
radiation  flux  away  from  a  planar  source.  Thus  the  flux-limiter  of  eqn.(12) 
should  be  correct  near  the  surface  of  a  cylindrical  source.  Using  eqn.(12) 
in  eqn.(7)  gives  a  diffusion  equation  for  the  radiation  energy  density; 


3E 

at 


3(c  +0  W2|  VE|/E 

3.  b 


VE 


ccj  E 
a 


e  . 


(13) 


This  equation  is  differenced  over  the  distorted  mesh  of  Fig.  2  in  a  nine 

point  stencil,  i.e.,  the  difference  equation  connects  nine  mesh  cells.  Nine 

cells  are  coupled  together  instead  of  the  usual  five  for  a  two  dimensional 

partial  differential  equation  due  to  the  non-rectangular  mesh.  The 
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differencing  scheme  is  similar  to  that  discussed  by  Kershaw  and  Pert 

During  a  time  step  the  diffusion  coefficient  in  eqn.(13)  is  held  constant  and 

the  resulting  set  of  linear  algebraic  equations  for  E  is  solved  by  an 

incomplete  Crout  decomposition  and  a  conjugate  gradient  iteration.  This 

19 

matrix  solver  is  similar  to  one  described  by  Hain  ,  but  with  improvements  to 


speed  convergence  and  generalized  to  allow  for  a  nine  point  stencil. 


The  numerical  solutions  for  the  radiation  energy  density  (ERAD)and 

normalized  radiation  flux  in  the  z-direction  (FZ/c*ERAD)  are  shown  in  Figs. 4a 

and  *b.  respectively.  The  radiation  flux  is  calculated  according  to  the 

::es  r.p'ron  for  the  flux-limiter,  eqn.(12),  after  the  solution  for  E  has 

•  The  presented  results  represent  twenty  calls  to  the  diffusion 

e:  -  .  • "  a  time  step  of  10  ^  seconds  per  call.  The  steady  state  solution 

-i  "  :er  o y  .noting  minimal  changes  from  the  solution  with  fewer  calls. 

' • : .  a . . . .  it  initially  takes  about  12  iterations  within  the  matrix  solver 

-4 

per  -all  to  reach  an  error  of  10  in  the  Lj  norm.  The  number  of  iterations 
drops  to  several  once  the  steady  state  is  approached.  The  boundary 
conditions  are  reflective  along  the  symmetry  axis  and  Dirichlet  with  E  =  0 
outside  the  remaining  boundaries. 


C.  Ray  Tracing  Solutions. 


The  essence  of  the  ray  tracing  approach  is  the  approximate  solution  of 
eqn.(4)  over  a  set  of  rays  extending  from  each  observation  point.  For  any 
mesh  we  take  the  observation  points  to  be  the  center  of  the  cells.  In  order 
to  have  a  general  approach  applicable  to  an  arbitrary  distribution  of  sources 
and  opacities  there  are  two  major  differences  from  the  closed  form  solution 
used  in  section  II. A. 

First,  the  optical  depth  and  line  integral  in  eqn.(4)  are  not  calculated 
exactly  as  for  the  closed  form  solution.  Instead  each  ray  from  an 
observation  point  r  is  broken  into  a  finite  number  of  equal  length  segments 
As  over  which  the  emission  and  absorption  coefficients  will  be  taken  to  be 
constant.  If  the  set  { s ^ ,  i  =  1  NSRAY  demarks  the  endpoints  of  the 
segments,  with  s.  >  s.+^,  the  solution  of  the  transrer  equation  over  a  single 
segment  along  the  direction  Q  can  be  written  as 


\\ 


I(r,s 


-a  .  ...As 

.+1,Q,v)  «  I(r,s.,Q,v)  e 

ei+l/2  (,  ~aa,  i+1/2^8 

+  l1 ' 6 

3  j 1  +  1/2 

Integration  along  the  ray  from  the  boundary  of  the  computational  volume  (s^  = 
s^)  to  the  observation  point  at  r  is  accomplished  by  repeated  application  of 
eqn.(14)  till  I(r, sNSRAY,  S,  u)  is  reached. 

Second,  for  the  closed  form  solution  only  those  rays  going  through  the 
source  region  were  considered.  In  the  general  ray  tracing  approach,  a  set  of 
direction  angles  (Cj)  j  =  1.  NRAY,  is  chosen  to  cover  half  the  unit  sphere 
about  an  observation  point,  and  this  set  is  the  same  for  each  observation 
point.  Only  half  the  unit  sphere  need  be  covered  due  to  the  symmetry  about 
the  computational  plane.  In  reference  to  Fig.  1,  there  are  NCOSTH  equally 
spaced  values  from  1  to  -1  for  the  cosine  of  the  polar  angle  0,  and  NPHI 
equally  spaced  angles  from  0  to  n  for  the  azimuthal  angle  <J>.  Thus  the  number 
of  rays  along  which  eqn.(14)  is  solved  is  NSRAY  =  (NC0STH-2)*NPHI  +  2.  Vhen 
account  is  taken  of  the  symmetry,  there  are  effectively  (NCOSTH-2)*(2*NPHI-2) 
+  2  rays.  The  weighting  assigned  to  each  ray  corresponds  to  a  solid  angle 
about  the  ray,  and  is  the  same  for  each  ray  except  for  the  two  polar 
directions. 

What  remains  to  be  discussed  is  the  method  for  evaluating  the  emissivity 

e.  ,  and  absorption  o  .  . for  each  segment  As  in  eqn.(14).  This  is  the 
1+1/2  a, 1+1/2 

heart  of  the  above  ray  trace  scheme  and  so  an  efficient  but  approximate 
method  was  developed.  In  Lagrangian  codes,  the  emissivity  and  absorption 
coefficients  are  taken  to  be  constant  within  each  cell  of  the  distorted  mesh. 
For  the  first  part  of  the  method,  a  mapping  from  the  distorted  mesh  in  the 
computational  plane  to  a  rectangular  mesh  is  formulated.  The  dimensions  of 


)  •  (14) 
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the  rectangular  mesh  are  arbitrary  with  Ar  the  width  of  each  rectangular 
cell,  Az  the  height.  Let  (K,L)  denote  the  distorted  mesh  and  (I,J)  the 
rectangular  mesh.  Then  for  an  arbitrary  cell  (I,J)  the  mapping  [  K  = 
KMAP(I,J),  L  =  LMAP(I,J)  ]  gives  the  cell  (K,L)  within  which  the  center  of 
cell  (I,J)  lies.  This  mapping  is  calculated  once  for  a  distorted  mesh. 
Next,  the  entire  computational  volume  is  embedded  in  a  X'-Y'-Z'  Cartesian 
coordinate  system  with  the  plane  of  the  computational  mesh  lying  in  the  Y'-Z' 
plane.  Now  for  each  segment  As  of  a  ray  in  eqn.(14)  the  corresponding 
direction  cosines  of  this  ray  can  be  used  to  determine  the  midpoint 
=  (X',Y',Z')  of  the  segment  As.  The  plasma  variables  at  this  midpoint  are 
the  same  as  those  at  the  point  (R  =  /(X'  +Y'  ),Z  =  Z' )  in  the  computational 
plane.  Division  of  R  by  Ar  and  Z  by  Az  then  gives  the  cell  (I,J)  within 
which  the  point  (R,Z)  lies,  and  finally  by  the  above  mapping,  a  computational 
cell  (K,L).  Once  this  computational  cell  has  been  determined  the  emissivity 
and  absorption  coefficients  for  the  segment  As  are  known.  It  is  obvious  that 
although  the  point  (R,Z)  lies  in  cell  (I,J),  and  the  center  of  cell  (I,J) 
lies  in  cell  (K,L),  it  is  not  necessarily  true  that  the  point  (R,Z)  lies  in 
cell  (K,L).  Hence,  improper  values  for  the  coefficients  could  be  used.  For 
gentle  gradients  in  the  coefficients  this  problem  would  not  be  severe,  while 
for  steep  gradients  an  increase  in  the  dimensions  of  the  rectangular  mesh 
would  minimize  the  errors.  For  the  present  problem  with  steep  gradients  we 
used  a  40x40  rectangular  mesh.  Clearly  the  cost  of  even  doubling  the  mesh  is 
minimal  in  storage  and  computation  time. 

The  results  of  the  ray  tracing  scheme  for  the  radiation  energy  density 
(ERAD)  and  normalized  radiation  flux  in  the  z-direction  (FZ/c*ERAD)  are  shown 
in  Figs. 5a  and  5b,  respectively.  The  parameters  NCOSTH  and  NPHI  for  the  set 
of  rays  are  shown  within  the  square  brackets  of  Fig. 5a.  The  radiation  flux 
is  calculated  from  the  defining  eqn.(8). 


1 1 


D.  Discussion. 


Let  us  first  compare  the  results  for  the  radiation  energy  density  in 
Figs. 3a,  4a,  and  5a.  The  quantity  in  parentheses  at  the  top  of  each  contour 
plot  is  the  total  radiation  energy  vithin  the  computational  volume.  The 
diffusion  result  is  high  by  -19%  and  the  ray  trace  result  is  lov  by  -7%, 
compared  to  the  closed  form  results.  This  tendency  is  reflected  in  the 
contour  plots  where  the  diffusive  solution  has  too  high  an  energy  densi  ty 
near  the  axis  and  in  the  outer  regions.  A  study  of  the  point  by  point  errors 
shows  that  the  fractional  error  of  the  diffusive  solution  increases  as  one 
moves  away  from  the  source  surface  to  over  +80%  near  R  =  1  cm,  while  the  ray 
trace  errors  are  somewhat  random  and  have  typical  extremes  of  ±20%.  The 
relative  smoothness  of  the  diffusion  result  compared  to  the  ray  trace  result 
is  due  to  the  fact  that  this  simple  test  problem  is  well  suited  for  a 
diffusion  approach:  the  entire  source  is  concentrated  in  a  region  and  the 
energy  density  naturally  decreases  away  from  this  region. 

A  comparison  of  the  normalized  radiation  flux  in  the  z-direction  from 
Figs.  3b,  4b,  and  5b  clearly  shows  the  inadequacy  of  a  diffusion  scheme  using 
a  flux-limiter  to  estimate  the  radiation  flux. 

To  show  the  effect  of  increasing  the  number  of  rays,  this  problem  was 
redone  using  the  ray  trace  approach  with  -10  times  more  rays,  but  with  the 
same  size  rectangular  (I,J)  mesh  and  number  of  segment  divisions  (NSRAY)  for 
each  ray  as  above.  The  result  for  the  radiation  energy  density  is  given  in 
Fig. 6,  which  is  nearly  identical  to  the  closed  form  solution. 
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E.  Further  Problems  and  Solutions. 


To  further  compare  the  diffusion  and  ray  trace  schemes  we  consider 
several  more  intricate  distributions  for  the  absorption  and  emission 
coefficients . 

For  the  next  test  problem  consider  the  same  configuration  as  the  above 
problem  but  add  a  purely  absorbing  torus  with  a  =  20.0  cm  ^  extending  from 
0.3  cm  to  0.4  cm  in  the  radial  direction  and  from  0.4  cm  to  0.6  cm  in  the  z- 
direction.  This  feature  should  create  a  shadow  immediately  behind  the 
absorbing  torus  away  from  the  source  core.  Ve  employ  a  rectangular  mesh  in 
the  computational  plane  for  this  and  the  following  problem  to  accomadate  the 
configurations.  The  solutions  from  the  diffusive  and  ray  trace  schemes  are 
shown  in  Figs. 7a  and  7b,  respectively.  Again  the  diffusion  scheme  gives  a 
higher  total  radiation  energy  in  the  computational  volume  than  the  ray  trace 
scheme.  More  significant  is  the  unphysical  character  of  the  energy  density 
contour  levels  for  the  diffusive  solution.  In  front  of  the  absorbing  torus, 
i.e.,  toward  the  source  region,  one  sees  in  Fig. 7a  that  the  40  and  35 
ergs/cm  -Hz  contours  are  attracted  toward  the  torus.  Behind  the  torus,  the 
expected  shadow  is  not  seen;  instead  one  finds  a  radiation  energy  densi  ty 

3 

between  5  and  10  ergs/cm  -Hz.  Note  for  the  ray  trace  solution  that  the 
contour  levels  have  a  reasonable  structure  in  front  of  the  torus,  and  the 
shadow  is  well  defined.  The  erroneous  features  of  the  diffusive  solution  are 
due  to  the  approximation  of  the  radiative  transfer  equation  as  a  diffusion 
equation.  In  a  diffusion  process  local  minimums  tend  to  fill  in.  For  this 
problem,  VE  at  the  edges  of  the  shadow  are  not  in  the  direction  of  the 
radiation  flux.  Hence  any  flux-limiting  diffusion  scheme  will  have 
difficulties  handling  shadows. 
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For  the  final  problem  let  an  emissivity  of  100  ergs/cm-  *secHz  be 
confined  to  two  square  tori  centered  about  the  symmetry  axis,  and  each 
extending  from  0.2  to  0.4  cm  in  radius.  In  the  z-direction,  the  lover  one 
extends  from  0.2  cm  to  0.4  cm,  and  the  upper  one  from  0.6  cm  to  0.8  cm. 
Outside  of  these  tori  the  emissivity  is  zero.  Let  the  absorption  coefficient 
be  o  =  1  cm  *  throughout  the  volume,  including  the  emitting  tori.  This 

d 

configuration  could  model  two  hot  spots  formed  during  a  collapsing  Z-pinch. 
The  radiation  energy  density  from  the  numerical  results  of  the  diffusion  and 
ray  trace  schemes  are  shown  in  Figs. 8a  and  8b,  respectively.  As  is  typical, 
the  diffusion  result  gives  ~50X  more  radiation  energy  integrated  throughout 


III.  CONCLUSIONS. 


We  have  presented  a  comparison  of  two  different  numerical  methods  for 
solving  the  radiative  transfer  equation  in  axisymmetric  geometries.  The 
solutions  for  these  tvo  schemes  were  compared  with  a  closed  form  solution  of 
a  simple,  steady  state  test  problem,  and  were  further  compared  against  each 
other  for  several  additional  problems  of  more  complex  spatial  configurations. 

There  are  several  potential  advantages  of  the  diffusive  approach  over 
the  ray  trace  approach.  First,  since  these  schemes  are  intended  for  use  in 
2-D,  cylindrical,  radiation-hydrodynamic  codes,  the  time  step  plays  a  role. 
For  the  implicit  diffusive  approach,  the  time  step  of  eqn.(13)  is  simply  that 
of  the  hydrodynamic  part  of  the  code.  The  present  ray  trace  scheme,  however, 
does  not  take  account  of  the  photon  travel  time.  For  laboratory  plasma 
experiments  as  mentioned  in  the  Introduction,  the  spatial  size  is  small 
enough  that  a  streaming  photon  will  leave  the  region  of  interest.  Under 
these  conditions  the  radiation  energy  density  is  small  compared  to  the  plasma 
thermal  energy  density,  except  possibly  in  optically  thick  regions.  The  time 
of  flight  problem  is  potentially  more  serious  in  many  as trophysical  problems 
due  to  the  large  length  scales.  In  such  cases  a  diffusion  approach  could 
follow  a  travelling  radiation  front. 

Second,  if  scattering  is  an  important  process,  the  diffusion  eqns.(7) 
and  (11)  can  readily  account  for  it.  This  is  clearly  not  so  for  the  ray 
trace  scheme  where  an  iteration  on  the  radiation  energy  density  would  need  to 
be  employed. 

The  diffusive  approach  can  also  accomodate  general  boundary  conditions, 
such  as  periodic  in  one  direction  or  reflective  at  a  boundary.  The  ray  trace 
scheme  is  limited  in  that  the  incident  intensity  upon  a  boundary  of  the 
computational  volume  must  be  known  a  priori .  This  limitation  is  negated  if 
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the  experimental  region  of  interest  can  be  enclosed  in  the  computational 
volume  without  reflective  boundaries. 

Fourth,  a  diffusive  approach  provides  smoother  solutions  than  a  ray 
trace  one  by  its  very  nature:  local  extrema  are  reduced  and  local  minima  are 
filled  in. 

However,  there  are  major  advantages  of  the  ray  trace  approach  over  the 
diffusive  one.  First,  as  the  test  problems  show,  the  present  diffusive 
scheme  can  overestimate  the  radiation  energy  density,  particularly  in 
optically  thin,  streaming  regions.  This  reflects  the  fact  that  in  such 
regions  the  radiation  diffusion  velocity  of  eqn.(12)  is  always  small  compared 
to  the  speed  of  light,  except  in  regions  of  steep  radiation  energy  gradients. 
Possible  corrections  for  this  problem  have  been  referenced  in  section  II. B. , 
but  they  all  have  fundamental  problems  in  purely  absorbing  regions.  The 
consequence  of  this  problem  for  radiation-hydrodynamic  codes  would  be  an 
overestimate  of  the  plasma  heating  rate  due  to  absorption.  Furthermore, 
shadows  caused  by  absorbing  spots  are  not  at  all  treated  correctly,  as 
evidenced  by  comparing  Figs. 7a  and  7b. 

Second,  unlike  a  diffusive  scheme,  a  general  ray  trace  scheme  as 

described  in  section  II. C.  can  easily  be  adjusted  to  check  for  accuracy  by 

increasing  the  number  of  rays.  On  the  other  hand,  the  number  of  rays  can  be 

reduced  to  speed  computation.  There  is  no  such  adaptibility  in  a  diffusive 

approach.  Ue  note  that  the  present  code  can  be  reduced  to  the  simpler  ray 
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trace  scheme  of  Colorabant  and  Vinsor  ’  by  limiting  the  number  of  rays  in  the 
azimuthal  direction  to  two  (NPHI=2). 

Finally,  the  present  diffusive  scheme  is  roughly  5  to  10  times  faster 
than  the  ray  trace  scheme  for  -65  rays  and  a  single  frequency  over  a  11x10 
computational  mesh.  However,  if  a  multi-frequency  calculation  is  required 
for  the  radiation-hydrodynamics  code,  then  the  ray  trace  scheme  would  be 


faster  for  more  than  -10  frequencies.  This  is  due  to  the  fact  that  in  the 
ray  trace  code  little  extra  computation  is  needed  for  more  than  one 
frequency:  in  integrating  the  transfer  equation  along  a  ray,  a  simple 
vectorizable  "do  loop"  over  the  frequencies  at  each  segment  As  suffices.  On 
the  other  hand,  for  any  diffusion  approach  the  diffusion  equation  must  be 
solved  over  the  whole  mesh  for  each  frequency. 
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Fig.  3  —  The  closed  form  solutions  to  the  simple  test  problem  of  a  finite  height  emitting  core,  (a)  the 
radiation  energy  density  in  ergs/cm1  Hz,  (b)  the  non-dimensional,  normalized  radiation  flux  in  the  z- 
direction 
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l  ig.  3  (Corn'd)  —  The  closed  form  solutions  to  the  simple  test  problem  of  a  finite  height  emitting  core; 
(a)  the  radiation  energy  density  in  ergs/cm’  •  Hz,  <b)  the  non-dimensional,  normalized  radiation  flux  in 
the  /-direction 
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Fig.  5  (Corn'd)  —  The  ray  trace  solutions  to  the  simple  test  problem;  (a)  the  radiation  energy  density 
in  ergs/cm3  •  Hz,  (b)  the  non-dimensional,  normalized  radiation  flux  in  the  z-direction. 
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Fig.  7  —  The  solutions  for  the  radiation  energy  density  of  the  second  test  problem  with  a  purely 
absorbing  torus  in  front  of  the  emitting  core;  (a)  the  diffusive  solution,  (b)  the  ray  trace  solution 
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—  The  solutions  for  the  radiation  energy  density  of  the  third  test  problem  with  two  emitting 
tori;  (a)  the  diffusive  solution,  (b)  the  ray  trace  solution. 
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